      subroutine nplotter_ttbar(p,wt,wt2,switch)
c--- Variable passed in to this routine:
c
c---      p:  4-momenta of particles in the format p(i,4)
c---          with the particles numbered according to the input file
c---          and components labelled by (px,py,pz,E)
c
c---     wt:  weight of this event
c
c---    wt2:  weight^2 of this event
c
c--- switch:  an integer equal to 0 or 1, depending on the type of event
c---                0  --> lowest order, virtual or real radiation
c---                1  --> counterterm for real radiation
      implicit none
      include 'vegas_common.f'
      include 'constants.f'
      include 'histo.f'
      include 'masses.f'
      include 'plabel.f'
      include 'process.f'
      double precision p(mxpart,4),wt,wt2,yrap,pt,ptjet(mxpart),
     & pttwo,ptthree,yrapthree,yrapsix,pt345678,y345678,m345678,mll,
     & ptl,yl,mwtest(mxpart),mt1,mt2,mwhad,pthad,yhad,ptttb,yttb,
     & ptfour,yrapfour,yrapseven,mt3,mt4,min14,min23,min24,dot,tiny,
     & ptsix,ptt,pttb,yt,ytb,mttb,mwp,mwm,mlb,ktdistbg,ktdistbbarg,
     & MpT,ht,sumE,sumX,sumY,sumZ,pxsum,pysum,pzsum,esum,sttj,rhottj,
     & HTvis, Meff, mET,mETx,mETy
      double precision mlbb,mlbbb,mlbmin,psquare,psquaredot,dphi
      double precision energy,piti,mtbest
      double precision plep(4),ptop(4),pleprest(4),ptoprest(4),mj1j2
      double precision pb(4),pbbar(4),sigmatot
      integer switch,n,nplotmax,j,
     & jetindex(mxpart),iorder(mxpart),ijet
      integer mypart,parta,partb,massij,mynpart
      character*4 tag
      logical first,creatent,dswhisto,dilepton,failed
      common/outputflags/creatent,dswhisto
      common/nplotmax/nplotmax
      parameter(tiny=1d-8)
      data first/.true./
      save first

************************************************************************
*                                                                      *
*     INITIAL BOOKKEEPING                                              *
*                                                                      *
************************************************************************

c--- determine whether or not process represents dilepton channel
      if ((case .eq. 'tt_bbl') .or. (case .eq. 'tt_ldk')) then
        dilepton=.true.
      else
        dilepton=.false.
      endif
      
      if (first) then
c--- Initialize histograms, without computing any quantities; instead
c--- set them to dummy values
        tag='book'
	mj1j2=0d0
        write(*,*) '**********************************'
        write(*,*) 'This is the first point, tag=',tag
        sigmatot=0d0
        write(*,*) 'sigmatot=',sigmatot
        write(*,*) '**********************************' 
        goto 99
      else
c--- Add event in histograms
        tag='plot'
      endif

************************************************************************
*                                                                      *
*     DEFINITIONS OF QUANTITIES TO PLOT                                *
*                                                                      *
************************************************************************

       ijet=0
       do j=2,9
         if ((plabel(j) .eq. 'pp') .or. (plabel(j) .eq. 'bq')
     &   .or.(plabel(j) .eq. 'ba')) then
           if (p(j,4) .gt. tiny) then
	     ijet=ijet+1
	     jetindex(ijet)=j
	     ptjet(ijet)=pt(j,p)
	   endif
         endif
       enddo   
       call arraysort(ijet,ptjet,iorder)

c--- Two hardest jets are now indexed by
c--- jetindex(iorder(1)) and jetindex(iorder(2))
       mj1j2=((p(jetindex(iorder(1)),4)+p(jetindex(iorder(2)),4))**2
     &       -(p(jetindex(iorder(1)),1)+p(jetindex(iorder(2)),1))**2
     &       -(p(jetindex(iorder(1)),2)+p(jetindex(iorder(2)),2))**2
     &       -(p(jetindex(iorder(1)),3)+p(jetindex(iorder(2)),3))**2)
       mj1j2=dsqrt(max(mj1j2,0d0))
       
************************************************************************
*                                                                      *
*     FILL HISTOGRAMS                                                  *
*                                                                      *
************************************************************************

c--- Call histogram routines
   99 continue

c--- Book and fill ntuple if that option is set, remembering to divide
c--- by # of iterations now that is handled at end for regular histograms
      if (creatent .eqv. .true.) then
        call bookfill(tag,p,wt/dfloat(itmx))  
      endif

c--- "n" will count the number of histograms
      n=1              

c--- Syntax of "bookplot" routine is:
c
c---   call bookplot(n,tag,titlex,var,wt,wt2,xmin,xmax,dx,llplot)
c
c---        n:  internal number of histogram
c---      tag:  "book" to initialize histogram, "plot" to fill
c---   titlex:  title of histogram
c---      var:  value of quantity being plotted
c---       wt:  weight of this event (passed in)
c---      wt2:  weight of this event (passed in)
c---     xmin:  lowest value to bin
c---     xmax:  highest value to bin
c---       dx:  bin width
c---   llplot:  equal to "lin"/"log" for linear/log scale

c--- This section tailored for comparison with CMS PAS TOP-11-013

      if (dilepton) then
c--- dilepton=specific plots

      sigmatot=sigmatot+wt
c      write(*,*) 'weight=',wt
c      write(*,*) 'sigmatot=',sigmatot
c
c*********************************
c ********* all the p**2 *********
c*********************************     

      psquare=dsqrt(
     &p(5,4)**2 - p(5,1)**2 - p(5,2)**2 - p(5,3)**2)
      call bookplot(n,tag,'sqrt(p^2(5))',psquare,wt,wt2,
     & 0d0,40d0,1d0,'lin')
      n=n+1

      psquaredot=dsqrt(max(dot(p,5,5),0d0))
      call bookplot(n,tag,'m(5)',psquaredot,wt,wt2,
     & 0d0,40d0,1d0,'lin')
      n=n+1

c*********************************
c ********* E=dsqrt(m^2+pT^2)*cosh[y] *********
c*********************************

c      piti=dsqrt( p(5,1)**2 + p(5,2)**2)
c      psquaredot=dsqrt(max(dot(p,5,5),0d0))
c      energy=dsqrt(piti**2+psquaredot**2)*cosh(yrap(5,p))
c      call bookplot(n,tag,'Et*cosh(y)_0_100',energy,wt,wt2,
c     & 0d0,100d0,1d0,'lin')
c      n=n+1

c-- *********************** 
      mypart=5
c-- ***********************      
      piti=dsqrt( p(mypart,1)**2 + p(mypart,2)**2)
      psquaredot=dsqrt(max(dot(p,mypart,mypart),0d0))
      energy=dsqrt(piti**2+psquaredot**2)*cosh(yrap(mypart,p))
      call bookplot(n,tag,'Et(5)*cosh(y)_0_100_mypart',energy,wt,wt2,
     & 0d0,100d0,1d0,'lin')
      n=n+1

      piti=dsqrt( p(mypart,1)**2 + p(mypart,2)**2)
      psquaredot=dsqrt(max(dot(p,mypart,mypart),0d0))
      energy=dsqrt(piti**2+psquaredot**2)*cosh(yrap(mypart,p))
      call bookplot(n,tag,'Et(5)*cosh(y)_100_200',energy,wt,wt2,
     & 100d0,200d0,1d0,'lin')
      n=n+1

      piti=dsqrt( p(mypart,1)**2 + p(mypart,2)**2)
      psquaredot=dsqrt(max(dot(p,mypart,mypart),0d0))
      energy=dsqrt(piti**2+psquaredot**2)*cosh(yrap(mypart,p))
      call bookplot(n,tag,'Et(5)*cosh(y)_200_300',energy,wt,wt2,
     & 200d0,300d0,1d0,'lin')
      n=n+1

      piti=dsqrt( p(mypart,1)**2 + p(mypart,2)**2)
      psquaredot=dsqrt(max(dot(p,mypart,mypart),0d0))
      energy=dsqrt(piti**2+psquaredot**2)*cosh(yrap(mypart,p))
      call bookplot(n,tag,'Et(5)*cosh(y)_300_400',energy,wt,wt2,
     & 300d0,400d0,1d0,'lin')
      n=n+1

c      piti=dsqrt( p(mypart,1)**2 + p(mypart,2)**2)
c      psquaredot=dsqrt(max(dot(p,mypart,mypart),0d0))
c      energy=dsqrt(piti**2+psquaredot**2)*cosh(yrap(mypart,p))
c      call bookplot(n,tag,'Et(6)*cosh(y)_0_400_4',energy,wt,wt2,
c     & 0d0,400d0,4d0,'lin')
c      n=n+1

c-- ***********************
      mypart=6
c-- ***********************
      piti=dsqrt( p(mypart,1)**2 + p(mypart,2)**2)
      psquaredot=dsqrt(max(dot(p,mypart,mypart),0d0))
      energy=dsqrt(piti**2+psquaredot**2)*cosh(yrap(mypart,p))
      call bookplot(n,tag,'Et(6)*cosh(y)_0_100_1',energy,wt,wt2,
     & 0d0,100d0,1d0,'lin')
      n=n+1

      piti=dsqrt( p(mypart,1)**2 + p(mypart,2)**2)
      psquaredot=dsqrt(max(dot(p,mypart,mypart),0d0))
      energy=dsqrt(piti**2+psquaredot**2)*cosh(yrap(mypart,p))
      call bookplot(n,tag,'Et(6)*cosh(y)_100_200_1',energy,wt,wt2,
     & 100d0,200d0,1d0,'lin')
      n=n+1

      piti=dsqrt( p(mypart,1)**2 + p(mypart,2)**2)
      psquaredot=dsqrt(max(dot(p,mypart,mypart),0d0))
      energy=dsqrt(piti**2+psquaredot**2)*cosh(yrap(mypart,p))
      call bookplot(n,tag,'Et(6)*cosh(y)_200_300_1',energy,wt,wt2,
     & 200d0,300d0,1d0,'lin')
      n=n+1

      piti=dsqrt( p(mypart,1)**2 + p(mypart,2)**2)
      psquaredot=dsqrt(max(dot(p,mypart,mypart),0d0))
      energy=dsqrt(piti**2+psquaredot**2)*cosh(yrap(mypart,p))
      call bookplot(n,tag,'Et(6)*cosh(y)_300_400_1',energy,wt,wt2,
     & 300d0,400d0,1d0,'lin')
      n=n+1

c      piti=dsqrt( p(mypart,1)**2 + p(mypart,2)**2)
c      psquaredot=dsqrt(max(dot(p,mypart,mypart),0d0))
c      energy=dsqrt(piti**2+psquaredot**2)*cosh(yrap(mypart,p))
c      call bookplot(n,tag,'Et(6)*cosh(y)_300_400_1',energy,wt,wt2,
c     & 0d0,400d0,1d0,'lin')
c      n=n+1

c-- **********************************************
c-- ****** mtop translation (mjet-wise) **********
c-- **********************************************

c      call bookplot(n,tag,'mb',mb,wt,wt2,
c     & 0d0,30d0,1d0,'lin')
c      n=n+1
c      call bookplot(n,tag,'mW',wmass,wt,wt2,
c     & 50d0,100d0,1d0,'lin')
c      n=n+1
c      write(*,*) 'mb=',mb
c      write(*,*) 'mW=',wmass
c  mW=   80.397999999999996
c  mb=   4.7500000000000000

c--- mtbest
      parta=5
      mtbest=p(parta,4)+dsqrt(p(parta,4)**2 - mb**2 + wmass**2)
      call bookplot(n,tag,'mtbest',mtbest,wt,wt2,
     & 50d0,500d0,5d0,'lin')
      parta=6
      mtbest=p(parta,4)+dsqrt(p(parta,4)**2 - mb**2 + wmass**2)
      call bookplot(n,tag,'mtbest',mtbest,wt,wt2,
     & 50d0,500d0,5d0,'lin')
      n=n+1

      parta=5
      mtbest=p(parta,4)+dsqrt(p(parta,4)**2 - mb**2 + wmass**2)
      call bookplot(n,tag,'mtbest_50_150_1',mtbest,wt,wt2,
     & 50d0,150d0,1d0,'lin')
      parta=6
      mtbest=p(parta,4)+dsqrt(p(parta,4)**2 - mb**2 + wmass**2)
      call bookplot(n,tag,'mtbest_50_150_1',mtbest,wt,wt2,
     & 50d0,150d0,1d0,'lin')
      n=n+1

      parta=5
      mtbest=p(parta,4)+dsqrt(p(parta,4)**2 - mb**2 + wmass**2)
      call bookplot(n,tag,'mtbest_150_250_1',mtbest,wt,wt2,
     & 150d0,250d0,1d0,'lin')
      parta=6
      mtbest=p(parta,4)+dsqrt(p(parta,4)**2 - mb**2 + wmass**2)
      call bookplot(n,tag,'mtbest_150_250_1',mtbest,wt,wt2,
     & 150d0,250d0,1d0,'lin')
      n=n+1

      parta=5
      mtbest=p(parta,4)+dsqrt(p(parta,4)**2 - mb**2 + wmass**2)
      call bookplot(n,tag,'mtbest_250_350_1',mtbest,wt,wt2,
     & 250d0,350d0,1d0,'lin')
      parta=6
      mtbest=p(parta,4)+dsqrt(p(parta,4)**2 - mb**2 + wmass**2)
      call bookplot(n,tag,'mtbest_250_350_1',mtbest,wt,wt2,
     & 250d0,350d0,1d0,'lin')
      n=n+1

      parta=5
      mtbest=p(parta,4)+dsqrt(p(parta,4)**2 - mb**2 + wmass**2)
      call bookplot(n,tag,'mtbest_350_450_1',mtbest,wt,wt2,
     & 350d0,450d0,1d0,'lin')
      parta=6
      mtbest=p(parta,4)+dsqrt(p(parta,4)**2 - mb**2 + wmass**2)
      call bookplot(n,tag,'mtbest_350_450_1',mtbest,wt,wt2,
     & 350d0,450d0,1d0,'lin')
      n=n+1

c--- mtbestfromWminus
      parta=7
      partb=8
      mtbest=p(parta,4)+p(partb,4)+dsqrt( (p(parta,4)+p(partb,4))**2 
     &+ mb**2 - wmass**2)
      call bookplot(n,tag,'mtbestWminus',mtbest,wt,wt2,
     & 50d0,500d0,5d0,'lin')
      n=n+1
c--- mtbestfromWplus
      parta=3
      partb=4
      mtbest=p(parta,4)+p(partb,4)+dsqrt( (p(parta,4)+p(partb,4))**2 
     &+ mb**2 - wmass**2)
      call bookplot(n,tag,'mtbestWplus',mtbest,wt,wt2,
     & 50d0,500d0,5d0,'lin')
      n=n+1

c--- mtbestfromWminus
      parta=7
      partb=8
      mtbest=p(parta,4)+p(partb,4)+dsqrt( (p(parta,4)+p(partb,4))**2
     &+ mb**2 - wmass**2)
      call bookplot(n,tag,'mtbestWminus_50_150_1',mtbest,wt,wt2,
     & 50d0,150d0,1d0,'lin')
      n=n+1
c--- mtbestfromWplus
      parta=3
      partb=4
      mtbest=p(parta,4)+p(partb,4)+dsqrt( (p(parta,4)+p(partb,4))**2
     &+ mb**2 - wmass**2)
      call bookplot(n,tag,'mtbestWplus_50_150_1',mtbest,wt,wt2,
     & 50d0,150d0,1d0,'lin')
      n=n+1

c--- mtbestfromWminus
      parta=7
      partb=8
      mtbest=p(parta,4)+p(partb,4)+dsqrt( (p(parta,4)+p(partb,4))**2
     &+ mb**2 - wmass**2)
      call bookplot(n,tag,'mtbestWminus_150_250_1',mtbest,wt,wt2,
     & 150d0,250d0,1d0,'lin')
      n=n+1
c--- mtbestfromWplus
      parta=3
      partb=4
      mtbest=p(parta,4)+p(partb,4)+dsqrt( (p(parta,4)+p(partb,4))**2
     &+ mb**2 - wmass**2)
      call bookplot(n,tag,'mtbestWplus_150_250_1',mtbest,wt,wt2,
     & 150d0,250d0,1d0,'lin')
      n=n+1

c--- mtbestfromWminus
      parta=7
      partb=8
      mtbest=p(parta,4)+p(partb,4)+dsqrt( (p(parta,4)+p(partb,4))**2
     &+ mb**2 - wmass**2)
      call bookplot(n,tag,'mtbestWminus_250_350_1',mtbest,wt,wt2,
     & 250d0,350d0,1d0,'lin')
      n=n+1
c--- mtbestfromWplus
      parta=3
      partb=4
      mtbest=p(parta,4)+p(partb,4)+dsqrt( (p(parta,4)+p(partb,4))**2
     &+ mb**2 - wmass**2)
      call bookplot(n,tag,'mtbestWplus_250_350_1',mtbest,wt,wt2,
     & 250d0,350d0,1d0,'lin')
      n=n+1

c--- mtbestfromWminus
      parta=7
      partb=8
      mtbest=p(parta,4)+p(partb,4)+dsqrt( (p(parta,4)+p(partb,4))**2
     &+ mb**2 - wmass**2)
      call bookplot(n,tag,'mtbestWminus_350_450_1',mtbest,wt,wt2,
     & 350d0,450d0,1d0,'lin')
      n=n+1
c--- mtbestfromWplus
      parta=3
      partb=4
      mtbest=p(parta,4)+p(partb,4)+dsqrt( (p(parta,4)+p(partb,4))**2
     &+ mb**2 - wmass**2)
      call bookplot(n,tag,'mtbestWplus_350_450_1',mtbest,wt,wt2,
     & 350d0,450d0,1d0,'lin')
      n=n+1


c-- **********************
c*********************************
c********** FM VARIABLES *********
c*********************************                   
c--- pt(l+ and l-)
c--- FM1
      call bookplot(n,tag,'pt(l+_and_l-)',pt(4,p),wt,wt2,
     & 0d0,400d0,10d0,'lin')
      call bookplot(n,tag,'pt(l+_and_l-)',pt(7,p),wt,wt2,
     & 0d0,400d0,10d0,'lin')
      n=n+1

c--- eta(l+ and l-)
c      call bookplot(n,tag,'eta(l+ and l-)',yrap(4,p),wt,wt2,
c     & -2.4d0,2.4d0,0.2d0,'lin')
c      call bookplot(n,tag,'eta(l+ and l-)',yrap(7,p),wt,wt2,
c     & -2.4d0,2.4d0,0.2d0,'lin')
c      n=n+1

c--- pt(l+,l-)
c--- FM2
      call bookplot(n,tag,'pt(l+,l-)',pttwo(4,7,p),wt,wt2,
     & 0d0,400d0,10d0,'lin')
      n=n+1

c--- m(l+,l-)
c--- FM3
      mll=dsqrt(max(2d0*dot(p,4,7),0d0))
      call bookplot(n,tag,'m(l+,l-)',mll,wt,wt2,
     & 0d0,400d0,10d0,'lin')
      n=n+1

c--- E(l-,l+)
c--- FM4
      call bookplot(n,tag,'E(l-,l+)',p(4,4)+p(7,4),wt,wt2,
     & 0d0,250d0,5d0,'lin')
      n=n+1

c--- pt(l+)+pt(l-)
c--- FM5
      call bookplot(n,tag,'pt(l+)+pt(l-)',pt(4,p)+pt(7,p),wt,wt2,
     & 0d0,400d0,10d0,'lin')
      n=n+1

c*********************************
c *********** ENERGIES ***********
c*********************************
c--- Brock proposal
c--- E(e- and e+)
      call bookplot(n,tag,'E(e-_and_e+)0_1000_20',p(4,4),wt,wt2,
     & 0d0,1000d0,20d0,'lin')
      call bookplot(n,tag,'E(e-_and_e+)0_1000_20',p(7,4),wt,wt2,
     & 0d0,1000d0,20d0,'lin')
      n=n+1

c--- E(e- and e+)
      call bookplot(n,tag,'E(e-_and_e+)0_400_4',p(4,4),wt,wt2,
     & 0d0,400d0,4d0,'lin')
      call bookplot(n,tag,'E(e-_and_e+)0_400_4',p(7,4),wt,wt2,
     & 0d0,400d0,4d0,'lin')
      n=n+1

c--- E(e- and e+)
      call bookplot(n,tag,'E(e-_and_e+)0_200_2',p(4,4),wt,wt2,
     & 0d0,200d0,2d0,'lin')
      call bookplot(n,tag,'E(e-_and_e+)0_200_2',p(7,4),wt,wt2,
     & 0d0,200d0,2d0,'lin')
      n=n+1

c ******************************
c ******  OUR OBSERVABLE *******
c ******************************
c--- E(b and bbar)
      call bookplot(n,tag,'E(b_and_bbar)',p(5,4),wt,wt2,
     & 0d0,400d0,4d0,'lin')
      call bookplot(n,tag,'E(b_and_bbar)',p(6,4),wt,wt2,
     & 0d0,400d0,4d0,'lin')
      n=n+1

c--- E(b and bbar)0_200_2
      call bookplot(n,tag,'E(b_and_bbar)0_200_2',p(5,4),wt,wt2,
     & 0d0,200d0,2d0,'lin')
      call bookplot(n,tag,'E(b_and_bbar)0_200_2',p(6,4),wt,wt2,
     & 0d0,200d0,2d0,'lin')
      n=n+1

c--- E(b and bbar)
      call bookplot(n,tag,'E(b_and_bbar)30_150_4',p(5,4),wt,wt2,
     & 30d0,150d0,4d0,'lin')
      call bookplot(n,tag,'E(b_and_bbar)30_150_4',p(6,4),wt,wt2,
     & 30d0,150d0,4d0,'lin')
      n=n+1

c--- E(b and bbar)
      call bookplot(n,tag,'E(b_and_bbar)1_101_1',p(5,4),wt,wt2,
     & 1d0,101d0,1d0,'lin')
      call bookplot(n,tag,'E(b_and_bbar)1_101_1',p(6,4),wt,wt2,
     & 1d0,101d0,1d0,'lin')
      n=n+1

c--- E(b and bbar)
      call bookplot(n,tag,'E(b_and_bbar)100_200_1',p(5,4),wt,wt2,
     & 100d0,200d0,1d0,'lin')
      call bookplot(n,tag,'E(b_and_bbar)100_200_1',p(6,4),wt,wt2,
     & 100d0,200d0,1d0,'lin')
      n=n+1

c--- E(b and bbar)
      call bookplot(n,tag,'E(b_and_bbar)200_300_1',p(5,4),wt,wt2,
     & 200d0,300d0,1d0,'lin')
      call bookplot(n,tag,'E(b_and_bbar)200_300_1',p(6,4),wt,wt2,
     & 200d0,300d0,1d0,'lin')
      n=n+1

c--- E(b and bbar)
      call bookplot(n,tag,'E(b_and_bbar)300_400_1',p(5,4),wt,wt2,
     & 300d0,400d0,1d0,'lin')
      call bookplot(n,tag,'E(b_and_bbar)300_400_1',p(6,4),wt,wt2,
     & 300d0,400d0,1d0,'lin')
      n=n+1


c--- **************************************
c--- ******** each single energy *********
c--- **************************************
c--- f(p1)+f(p2) -->
c---  t(-->nu(p3)+e^+(p4)+b(p5))
c---  t~(-->b~(p6)+e^-(p7)+nu~(p8))

c--- E(1)
c    call bookplot(n,tag,'E(1)',p(1,4),wt,wt2,
c     & 0d0,250d0,5d0,'lin')
c      n=n+1
c--- E(2)
c      call bookplot(n,tag,'E(2)',p(2,4),wt,wt2,
c     & 0d0,250d0,5d0,'lin')
c      n=n+1
c
c
c--- E(3)
      call bookplot(n,tag,'E(nu)',p(3,4),wt,wt2,
     & 0d0,250d0,5d0,'lin')
      n=n+1
c--- E(4)
      call bookplot(n,tag,'E(e+)',p(4,4),wt,wt2,
     & 0d0,250d0,5d0,'lin')
      n=n+1
      call bookplot(n,tag,'E(e+)_0_100_1',p(4,4),wt,wt2,
     & 0d0,100d0,1d0,'lin')
      n=n+1
      call bookplot(n,tag,'E(e+)_100_200_1',p(4,4),wt,wt2,
     & 100d0,200d0,1d0,'lin')
      n=n+1
      call bookplot(n,tag,'E(e+)_200_300_1',p(4,4),wt,wt2,
     & 200d0,300d0,1d0,'lin')
      n=n+1
      call bookplot(n,tag,'E(e+)_300_400_1',p(4,4),wt,wt2,
     & 300d0,400d0,1d0,'lin')
      n=n+1
c--- E(5)
      call bookplot(n,tag,'E(b)',p(5,4),wt,wt2,
     & 0d0,250d0,5d0,'lin')
      n=n+1
      call bookplot(n,tag,'E(b)_0_100_1',p(5,4),wt,wt2,
     & 0d0,100d0,1d0,'lin')
      n=n+1
      call bookplot(n,tag,'E(b)_100_200_1',p(5,4),wt,wt2,
     & 100d0,200d0,1d0,'lin')
      n=n+1
      call bookplot(n,tag,'E(b)_200_300_1',p(5,4),wt,wt2,
     & 200d0,300d0,1d0,'lin')
      n=n+1
      call bookplot(n,tag,'E(b)_300_400_1',p(5,4),wt,wt2,
     & 300d0,400d0,1d0,'lin')
      n=n+1
c--- E(6)
      call bookplot(n,tag,'E(bbar)',p(6,4),wt,wt2,
     & 0d0,250d0,5d0,'lin')
      n=n+1
      call bookplot(n,tag,'E(bbar)_0_100_1',p(6,4),wt,wt2,
     & 0d0,100d0,1d0,'lin')
      n=n+1
      call bookplot(n,tag,'E(bbar)_100_200_1',p(6,4),wt,wt2,
     & 100d0,200d0,1d0,'lin')
      n=n+1
      call bookplot(n,tag,'E(bbar)_200_300_1',p(6,4),wt,wt2,
     & 200d0,300d0,1d0,'lin')
      n=n+1
      call bookplot(n,tag,'E(bbar)_300_400_1',p(6,4),wt,wt2,
     & 300d0,400d0,1d0,'lin')
      n=n+1
c--- E(7)
      call bookplot(n,tag,'E(e-)',p(7,4),wt,wt2,
     & 0d0,250d0,5d0,'lin')
      n=n+1
      call bookplot(n,tag,'E(e-)_0_100_1',p(7,4),wt,wt2,
     & 0d0,100d0,1d0,'lin')
      n=n+1
      call bookplot(n,tag,'E(e-)_100_200_1',p(7,4),wt,wt2,
     & 100d0,200d0,1d0,'lin')
      n=n+1
      call bookplot(n,tag,'E(e-)_200_300_1',p(7,4),wt,wt2,
     & 200d0,300d0,1d0,'lin')
      n=n+1
      call bookplot(n,tag,'E(e-)_300_400_1',p(7,4),wt,wt2,
     & 300d0,400d0,1d0,'lin')
      n=n+1
c--- E(8)
      call bookplot(n,tag,'E(nubar)',p(8,4),wt,wt2,
     & 0d0,250d0,5d0,'lin')
      n=n+1
c--- E(9)
      call bookplot(n,tag,'E(9)',p(9,4),wt,wt2,
     & 0d0,250d0,5d0,'lin')
      n=n+1
      call bookplot(n,tag,'E(9)_0_100_1',p(9,4),wt,wt2,
     & 0d0,100d0,1d0,'lin')
      n=n+1
      call bookplot(n,tag,'E(9)_100_200_1',p(9,4),wt,wt2,
     & 100d0,200d0,1d0,'lin')
      n=n+1
      call bookplot(n,tag,'E(9)_200_300_1',p(9,4),wt,wt2,
     & 200d0,300d0,1d0,'lin')
      n=n+1
      call bookplot(n,tag,'E(9)_300_400_1',p(9,4),wt,wt2,
     & 300d0,400d0,1d0,'lin')
      n=n+1
      call bookplot(n,tag,'pt(9)_0_400_4',pt(9,p),wt,wt2,
     & 0d0,400d0,4d0,'lin')
      n=n+1
c      
c--- ***** dphi for kT distance ****
      parta=9
      partb=5
      dphi=acos( (p(parta,1)*p(partb,1) +
     & p(parta,2)*p(partb,2))/
     & (dsqrt(p(parta,1)**2 + p(parta,2)**2 )*
     &dsqrt(p(partb,1)**2 + p(partb,2)**2 )))
c---
      ktdistbg=min(pt(parta,p),pt(partb,p))*dsqrt( dphi**2
     &+ (yrap(parta,p)-yrap(partb,p))**2 )
      call bookplot(n,tag,'kT(9,5)',kTdistbg,wt,wt2,
     & 0d0,80d0,1d0,'lin')
      n=n+1
c--- ***** dphi for kT distance ****
      parta=9
      partb=6
      dphi=acos( (p(parta,1)*p(partb,1) +
     & p(parta,2)*p(partb,2))/
     & (dsqrt(p(parta,1)**2 + p(parta,2)**2 )*
     &dsqrt(p(partb,1)**2 + p(partb,2)**2 )))
c---
      ktdistbg=min(pt(parta,p),pt(partb,p))*dsqrt( dphi**2
     &+ (yrap(parta,p)-yrap(partb,p))**2 )
      call bookplot(n,tag,'kT(9,6)',kTdistbg,wt,wt2,
     & 0d0,80d0,1d0,'lin')
      n=n+1


c*******************************************
c************* BMS 1006.0910 ***************
c*******************************************

c--- FIG.7 of 1006.0910 
c--- E(b,bbar)
      call bookplot(n,tag,'E(b,bbar)',p(5,4)+p(6,4),wt,wt2,
     & 0d0,1500d0,20d0,'lin')
      n=n+1

c--- FIG.5 of 1006.0910 
c--- E(e-,e+)
      call bookplot(n,tag,'E(e-,e+)',p(4,4)+p(7,4),wt,wt2,
     & 0d0,1000d0,20d0,'lin')
      n=n+1


c--- FIG.3 of 1006.0910
c--- min ( m(l+,b), m(l+,bbar)) lin
      mlbb=dsqrt((p(4,4)+p(5,4))**2-(p(4,1)+p(5,1))**2
     &         -(p(4,2)+p(5,2))**2-(p(4,3)+p(5,3))**2)
      mlbbb=dsqrt((p(4,4)+p(6,4))**2-(p(4,1)+p(6,1))**2
     &         -(p(4,2)+p(6,2))**2-(p(4,3)+p(6,3))**2)

      mlbmin = -100d0

         if ( mlbb.le.mlbbb  ) then
           mlbmin=mlbb
         endif

         if ( mlbbb.lt.mlbb  ) then
           mlbmin=mlbbb
         endif
            
      call bookplot(n,tag,'min(mlb)',mlbmin,wt,wt2,
     & 0d0,170d0,5d0,'lin')
      n=n+1
       

c--- eta(l+ and l-)
c      call bookplot(n,tag,'eta(l+ and l-)',yrap(4,p),wt,wt2,
c     & -2.4d0,2.4d0,0.2d0,'lin')
c      call bookplot(n,tag,'eta(l+ and l-)',yrap(7,p),wt,wt2,
c     & -2.4d0,2.4d0,0.2d0,'lin')
c      n=n+1

c--- pt(l+,l-)
c      call bookplot(n,tag,'pt(l+,l-)',pttwo(4,7,p),wt,wt2,
c     & 0d0,400d0,10d0,'log')
c      n=n+1

c--- m(l+,l-)
c      mll=dsqrt(max(2d0*dot(p,4,7),0d0))
c      call bookplot(n,tag,'m(l+,l-)',mll,wt,wt2,
c     & 0d0,400d0,10d0,'log')
c      n=n+1
                  
      else

c--- lepton+jets-specific plots
      if (plabel(4) .eq. 'ea') then
        ptl=pt(4,p)
	yl=yrap(4,p)
      else
        ptl=pt(7,p)
	yl=yrap(7,p)
      endif
      
c--- pt(lepton)
c      call bookplot(n,tag,'pt(lepton)',ptl,wt,wt2,
c     & 0d0,200d0,10d0,'log')
c      n=n+1

c--- eta(lepton)
c      call bookplot(n,tag,'eta(lepton)',yl,wt,wt2,
c     & -2.2d0,2.2d0,0.2d0,'lin')
c      n=n+1

      endif

c--- these plots are common to dilepton and lepton+jets channels
      call topreconstruct(p,failed,
     & ptt,yt,pttb,ytb,ptttb,yttb,mttb,mwp,mwm)

      if ((failed) .and. (first .eqv. .false.)) then
c--- make sure to increment by the number of plots in the 'else' section
      n=n+15
      
      else
      
c--- pt(t)
      call bookplot(n,tag,'pt(t)',ptt,wt,wt2,
     & 0d0,400d0,10d0,'log')
      n=n+1
      
c--- pt(tbar)
      call bookplot(n,tag,'pt(tbar)',pttb,wt,wt2,
     & 0d0,400d0,10d0,'log')
      n=n+1

c--- pt(t and tbar)
      call bookplot(n,tag,'pt(t_and_tbar)',ptt,wt,wt2,
     & 0d0,400d0,10d0,'log')
      call bookplot(n,tag,'pt(t_and_tbar)',pttb,wt,wt2,
     & 0d0,400d0,10d0,'log')
      n=n+1

c--- y(t)
      call bookplot(n,tag,'y(t)',yt,wt,wt2,
     & -2.6d0,2.6d0,0.2d0,'lin')
      n=n+1

c--- y(tbar)
      call bookplot(n,tag,'y(tbar)',ytb,wt,wt2,
     & -2.6d0,2.6d0,0.2d0,'lin')
      n=n+1

c--- y(t and tbar)
      call bookplot(n,tag,'y(t_and_tbar)',yt,wt,wt2,
     & -2.6d0,2.6d0,0.2d0,'lin')
      call bookplot(n,tag,'y(t_and_tbar)',ytb,wt,wt2,
     & -2.6d0,2.6d0,0.2d0,'lin')
      n=n+1

c--- pt(t,tbar)
      call bookplot(n,tag,'pt(t,tbar)',ptttb,wt,wt2,
     & -0.01d0,399.99d0,10d0,'log')
      n=n+1

c-- ***********************
c--        s(tt+j)
c-- ***********************

      if (abs(p(9,4)) .gt. 1d-8) then
        mynpart=9
      else
        mynpart=8
      endif

      pxsum=0d0
      pysum=0d0
      pzsum=0d0
      esum=0d0
      do j=3,mynpart
      pxsum=pxsum+p(j,1)
      pysum=pysum+p(j,2)
      pzsum=pzsum+p(j,3)
      esum=esum+p(j,4)
      enddo

      sttj=dsqrt(esum*esum - pxsum*pxsum  
     & - pysum*pysum - pzsum*pzsum)

      rhottj=340.0d0/sttj

c--- rhottj(3,...,npart-2)
      call bookplot(n,tag,'rhottj',rhottj,wt,wt2,
     & 0.0d0,1.00d0,0.025d0,'lin')
      n=n+1

c--- sttj(3,...,npart-2)
      call bookplot(n,tag,'sttj_240_340',sttj,wt,wt2,
     & 240.0d0,340.00d0,1d0,'lin')
      n=n+1

c--- sttj(3,...,npart-2)
      call bookplot(n,tag,'sttj_340_2140',sttj,wt,wt2,
     & 340.0d0,2140.00d0,45d0,'lin')
      n=n+1

c--- sttj(3,...,npart-2)
      call bookplot(n,tag,'sttj_340_640',sttj,wt,wt2,
     & 340.0d0,640.00d0,7.5d0,'lin')
      n=n+1

c--- sttj(3,...,npart-2)
      call bookplot(n,tag,'sttj_640_940',sttj,wt,wt2,
     & 640.0d0,940.00d0,7.5d0,'lin')
      n=n+1

c--- sttj(3,...,npart-2)
      call bookplot(n,tag,'sttj_940_1540',sttj,wt,wt2,
     & 940.0d0,1540.00d0,15d0,'lin')
      n=n+1

c--- sttj(3,...,npart-2)
      call bookplot(n,tag,'sttj_1540_2140',sttj,wt,wt2,
     & 1540.0d0,2140.00d0,15d0,'lin')
      n=n+1

c-- ***********************
c--        HT
c-- ***********************

      if (abs(p(9,4)) .gt. 1d-8) then
        mynpart=9
      else
        mynpart=8
      endif

      ht=0d0
      do j=3,mynpart
      ht=ht+pt(j,p)
      enddo

c--- HT(3,...,npart-2)
      call bookplot(n,tag,'HT_0_600_15',ht,wt,wt2,
     & 0.0d0,600.00d0,15d0,'lin')
      n=n+1

c--- HT(3,...,npart-2)
      call bookplot(n,tag,'HT_600_1200_15',ht,wt,wt2,
     & 600.0d0,1200.00d0,15d0,'lin')
      n=n+1

c--- HT(3,...,npart-2)
      call bookplot(n,tag,'HT_0_1200_30',ht,wt,wt2,
     & 0.0d0,1200.00d0,30d0,'lin')
      n=n+1

c-- ***********************
c--        HTivis (emubb) 4,5,6,7
c--        mET   (nu nu ) 3,8
c--        Meff
c-- ***********************
      if (abs(p(9,4)) .gt. 1d-8) then
        mynpart=9
      else
        mynpart=8
      endif

      HTvis=0.0d0
      mET=0.0d0
      mETx=0.0d0
      mETy=0.0d0
      Meff=0.0d0
      do j=3,mynpart
        if ( ( j .eq. 3 ) .or. ( j .eq. 8  ) ) then
                 mETx=mETx+p(j,1)
                 mETy=mETy+p(j,2)
        else 
                 HTvis=HTvis+pt(j,p)
        endif
      enddo
      mET = dsqrt( mETx**2 + mETy**2 )
      Meff = mET + HTvis   

c--- HTvis
      call bookplot(n,tag,'HTvis_0_400_4',HTvis,wt,wt2,
     & 0.0d0,400.00d0,4d0,'lin')
      n=n+1
c--- HTvis
      call bookplot(n,tag,'HTvis_400_800_4',HTvis,wt,wt2,
     & 400.0d0,800.00d0,4d0,'lin')
      n=n+1
c--- HTvis
      call bookplot(n,tag,'HTvis_800_1200_4',HTvis,wt,wt2,
     & 800.0d0,1200.00d0,4d0,'lin')
      n=n+1
c--- HTvis
      call bookplot(n,tag,'HTvis_1200_2400_12',HTvis,wt,wt2,
     & 1200.0d0,2400.00d0,12d0,'lin')
      n=n+1
c--- mET
      call bookplot(n,tag,'mET_0_400_4',mET,wt,wt2,
     & 0.0d0,400.00d0,4d0,'lin')
      n=n+1
c--- mET
      call bookplot(n,tag,'mET_400_800_4',mET,wt,wt2,
     & 400.0d0,800.00d0,4d0,'lin')
      n=n+1
c--- mET
      call bookplot(n,tag,'mET_800_1200_4',mET,wt,wt2,
     & 800.0d0,1200.00d0,4d0,'lin')
      n=n+1
c--- mET
      call bookplot(n,tag,'mET_1200_2400_12',mET,wt,wt2,
     & 1200.0d0,2400.00d0,12d0,'lin')
      n=n+1

c--- Meff
      call bookplot(n,tag,'Meff_0_400_4',Meff,wt,wt2,
     & 0.0d0,400.00d0,4d0,'lin')
      n=n+1
c--- Meff
      call bookplot(n,tag,'Meff_400_800_4',Meff,wt,wt2,
     & 400.0d0,800.00d0,4d0,'lin')
      n=n+1
c--- Meff
      call bookplot(n,tag,'Meff_800_1200_4',Meff,wt,wt2,
     & 800.0d0,1200.00d0,4d0,'lin')
      n=n+1
c--- Meff
      call bookplot(n,tag,'Meff_1200_2400_12',Meff,wt,wt2,
     & 1200.0d0,2400.00d0,12d0,'lin')
      n=n+1


c-- ***********************
c--        MpT
c-- ***********************
      sumX=p(3,1) + p(4,1) + p(5,1)
      sumY=p(3,2) + p(4,2) + p(5,2)
      sumZ=p(3,3) + p(4,3) + p(5,3)     
      sumE=p(3,4) + p(4,4) + p(5,4)

      piti=dsqrt( sumX**2 + sumY**2 )
      MpT=dsqrt( sumE**2 - 
     & sumX**2 - sumY**2 - sumZ**2 + piti**2 )
c---     MpT=dsqrt( sumE**2 - sumZ**2 )

c--- MpT(3,4,5)
      call bookplot(n,tag,'MpT(3,4,5)_0_400_10',MpT,wt,wt2,
     & 0.0d0,400.00d0,10d0,'lin')
      n=n+1

c--- MpT(3,4,5)
      call bookplot(n,tag,'MpT(3,4,5)_400_800_10',MpT,wt,wt2,
     & 400.0d0,800.00d0,10d0,'lin')
      n=n+1

c--- MpT(3,4,5)
      call bookplot(n,tag,'MpT(3,4,5)_0_800_20',MpT,wt,wt2,
     & 0.0d0,800.00d0,20d0,'lin')
      n=n+1

c-- ***********************
c-- ***********************
c-- ***********************

c--- y(t,tbar)
      call bookplot(n,tag,'y(t,tbar)',yttb,wt,wt2,
     & -2.6d0,2.6d0,0.2d0,'lin')
      n=n+1

c--- m(t,tbar)
      call bookplot(n,tag,'m(t,tbar)log',mttb,wt,wt2,
     & 0d0,1400d0,25d0,'log')
      n=n+1

c--- m(t,tbar)
      call bookplot(n,tag,'m(t,tbar)lin',mttb,wt,wt2,
     & 0d0,1400d0,25d0,'lin')
      n=n+1

c--- m(W+ candidate)
      call bookplot(n,tag,'m(W+_candidate)lin',mwp,wt,wt2,
     & 0d0,200d0,5d0,'lin')
      n=n+1

c--- m(W- candidate)
      call bookplot(n,tag,'m(W-_candidate)lin',mwm,wt,wt2,
     & 0d0,200d0,5d0,'lin')
      n=n+1

c--- m(W+ candidate)
      call bookplot(n,tag,'m(W+_candidate)log',mwp,wt,wt2,
     & 0d0,200d0,5d0,'log')
      n=n+1

c--- m(W- candidate)
      call bookplot(n,tag,'m(W-_candidate)log',mwm,wt,wt2,
     & 0d0,200d0,5d0,'log')
      n=n+1
      endif
      
c--- This section tailored for comparison with CMS PAS TOP-11-013

c--- plots for comparison with arXiv:0907.3090
      if (dilepton) then

c--- pt(l+)
c      call bookplot(n,tag,'pt(l+)',pt(4,p),wt,wt2,
c     & 0d0,1000d0,50d0,'log')
c      n=n+1
c--- pt(l-)
c      call bookplot(n,tag,'pt(l-)',pt(7,p),wt,wt2,
c     & 0d0,1000d0,50d0,'log')
c      n=n+1

c--- eta(l+)
c      call bookplot(n,tag,'eta(l+)',yrap(4,p),wt,wt2,
c     & -4d0,4d0,0.2d0,'lin')
c      n=n+1
c--- eta(l-)
c      call bookplot(n,tag,'eta(l-)',yrap(7,p),wt,wt2,
c     & -4d0,4d0,0.2d0,'lin')
c      n=n+1
      endif

c--- m(l+,b)  log
      mlb=dsqrt((p(4,4)+p(5,4))**2-(p(4,1)+p(5,1))**2
     &         -(p(4,2)+p(5,2))**2-(p(4,3)+p(5,3))**2)
      call bookplot(n,tag,'m(l+,b)log',mlb,wt,wt2,
     & 0d0,200d0,5d0,'log')
      n=n+1
c--- m(l-,bb)   log
      mlb=dsqrt((p(6,4)+p(7,4))**2-(p(6,1)+p(7,1))**2
     &         -(p(6,2)+p(7,2))**2-(p(6,3)+p(7,3))**2)
      call bookplot(n,tag,'m(l-,bb)log',mlb,wt,wt2,
     & 0d0,200d0,5d0,'log')
      n=n+1
c--- m(l+,b)   lin
      mlb=dsqrt((p(4,4)+p(5,4))**2-(p(4,1)+p(5,1))**2
     &         -(p(4,2)+p(5,2))**2-(p(4,3)+p(5,3))**2)
      call bookplot(n,tag,'m(l+,b)lin',mlb,wt,wt2,
     & 0d0,200d0,5d0,'lin')
      n=n+1
c--- m(l-,bb)   lin
      mlb=dsqrt((p(6,4)+p(7,4))**2-(p(6,1)+p(7,1))**2
     &         -(p(6,2)+p(7,2))**2-(p(6,3)+p(7,3))**2)
      call bookplot(n,tag,'m(l-,bb)lin',mlb,wt,wt2,
     & 0d0,200d0,5d0,'lin')
      n=n+1
c--- plots for comparison with arXiv:0907.3090
      

c--- plots for computing ttbar asymmetry
  
      if ((failed) .and. (first .eqv. .false.)) then
c--- make sure to increment by the number of plots in the 'else' section
      n=n+3
      
      else
            
c--- inclusive case
c      call bookplot(n,tag,'y(t)-y(tbar)',yt-ytb,wt,wt2,
c     & -4d0,4d0,4d0,'lin')
c      n=n+1

c--- for bins of |y(t)-y(tbar)| of (0,0.5), (0.5,1), (1,1.5)
c      call bookplot(n,tag,'y(t)-y(tbar)',yt-ytb,wt,wt2,
c     & -1.5d0,1.5d0,0.5d0,'lin')
c      n=n+1

c--- for bins of |y(t)-y(tbar)| of (1.5,4.5)
c      call bookplot(n,tag,'y(t)-y(tbar)',yt-ytb,wt,wt2,
c     & -4.5d0,4.5d0,3d0,'lin')
c      n=n+1

      endif
      
c-- lepton asymmetry
      if (dilepton .eqv. .false.) then
        if (plabel(4) .eq. 'ea') then
c        call bookplot(n,tag,'Q(lep)*y(lep)',yrap(4,p),wt,wt2,
c     &   -4d0,4d0,4d0,'lin')
c        n=n+1
	endif
        if (plabel(7) .eq. 'el') then
c        call bookplot(n,tag,'Q(lep)*y(lep)',-yrap(7,p),wt,wt2,
c     &   -4d0,4d0,4d0,'lin')
c        n=n+1
	endif
      endif
      
c--- plots for computing ttbar asymmetry

c--- invariant mass of two leading jets       
c      call bookplot(n,tag,'m(j1,j2)',mj1j2,wt,wt2, 0d0,300d0,5d0,'lin')
c      n=n+1
      
c--- compute energy of particle 4 in rest frame of top
      do j=1,4
      plep(j)=p(4,j)
      ptop(j)=p(3,j)+p(4,j)+p(5,j)
      ptoprest(j)=0d0
      enddo
      ptoprest(4)=mt
c      call boostx(plep,ptop,ptoprest,pleprest)      
c      call bookplot(n,tag,'2*E(4)/mt (top rest frame)',
c     & 2d0*pleprest(4)/mt,wt,wt2,0d0,1d0,0.01d0,'lin')
c      n=n+1

c--- compute energy of particle 7 in rest frame of top
      do j=1,4
      plep(j)=p(7,j)
      ptop(j)=p(6,j)+p(7,j)+p(8,j)
      ptoprest(j)=0d0
      enddo
      ptoprest(4)=mt
c      call boostx(plep,ptop,ptoprest,pleprest)      
c      call bookplot(n,tag,'2*E(7)/mt (top rest frame)',
c     & 2d0*pleprest(4)/mt,wt,wt2,0d0,1d0,0.01d0,'lin')
c      n=n+1
       
c--- single-particle plots
      do j=3,8
        call genplot1(p,j,tag,wt,wt2,n)
      enddo
c--- two-particle plots
      call genplot2(p,3,4,tag,wt,wt2,n)
      call genplot2(p,7,8,tag,wt,wt2,n)
      call genplot2(p,5,6,tag,wt,wt2,n)
c--- three-particle plots
      call genplot3(p,3,4,5,tag,wt,wt2,n)
      call genplot3(p,6,7,8,tag,wt,wt2,n)

c--- additional plots that may be present at NLO       
      if (abs(p(9,4)) .gt. 1d-8) then
        call genplot1(p,9,tag,wt,wt2,n)
        call genplot3(p,7,8,9,tag,wt,wt2,n)
      else
        n=n+2
      endif

************************************************************************
*                                                                      *
*     FINAL BOOKKEEPING                                                *
*                                                                      *
************************************************************************

c--- We have over-counted the number of histograms by 1 at this point
      n=n-1

c--- Ensure the built-in maximum number of histograms is not exceeded    
      call checkmaxhisto(n)

c--- Set the maximum number of plots, on the first call
      if (first) then
        first=.false.
        nplotmax=n
      endif
      
      return
      end
